import numpy as np
from scipy.stats import norm


class BSM:
    """Black Scholes Model
    Example:

    K = 100;  r = 0.1 ;   T = 1 ;   sigma = 0.3 ;   S = np.arange(60,140,0.1)

    https://www.codearmo.com/python-tutorial/options-trading-greeks-black-scholes

    Attributes:
        S : current asset price
        K: strike price of the option
        r: risk free rate
        T : time until option expiration
        σ: annualized volatility of the asset's returns
        q: continuously compounded dividend ,0 to ignore dividend

    Returns:
        Call: Call=S0*N(d1)*(e^-qT)-N(d2)K*(e^-rT)
        Put: Put=N(-d2)*K*(e^-rT)-N(-d1)*S0*(e^-qT)
        d1: d1=[ln(S/K)+(r-q+omega^2/2)*T]/[omega*(T^1/2)]
        d2: d2=d1-omega*(T^1/2)
    Raises:
        cpt err: every time new BSM instance to avoid grade cpt pollute attribution value.
        use greeks: https://zhuanlan.zhihu.com/p/90290453


    """

    def __init__(self, S, K, T, r, opt_ty="call", sigma=0.5, q=0.0):
        self.S = S
        self.K = K
        self.T = T
        self.r = r
        self.opt_ty = opt_ty
        self.sigma = sigma
        self.q = q
        self.ds = 1e-4

    def N(self, x):
        return norm.cdf(x)

    def dN(self, x):
        return np.exp(-x * x / 2) / np.sqrt(2 * np.pi)

    @property
    def factor_q(self):
        return np.exp(-self.q * self.T)

    @property
    def factor_r(self):
        return np.exp(-self.r * self.T)

    @property
    def d1(self):
        up = np.log(self.S / self.K) + (self.r - self.q + self.sigma ** 2 / 2) * self.T
        return up / (self.sigma * np.sqrt(self.T))

    @property
    def d2(self):
        return self.d1 - self.sigma * np.sqrt(self.T)

    def call_value(self):
        return self.S * self.factor_q * self.N(self.d1) - \
               self.K * self.factor_r * self.N(self.d2)

    def put_value(self):
        return self.K * self.factor_r * self.N(-self.d2) - \
               self.S * self.factor_q * self.N(-self.d1)

    @property
    def price(self):
        if self.opt_ty == 'call':
            return self.call_value()
        if self.opt_ty == 'put':
            return self.put_value()
        else:
            raise ValueError('Unrecognized type')

    def implied_sigma(self, market_price, one_cent=0.01, max_try=1000):
        """隐含波动率  牛顿逼近 本质就是梯度下降,从self.sigma起步搜索 默认0.5
            Args:
                option_ty: C for Call ; P for Put
                market_price: market price
                sigma0: start point for sigma
                max_try: how many loops to try
        """
        for i in range(max_try):
            diff = market_price - self.price
            if abs(diff) < one_cent:
                return self.sigma
            self.sigma += diff / self.vega
        return self.sigma

    @property
    def delta(self):
        """基础价格对期权价格求导
        """
        nd1 = self.N(self.d1)
        if self.opt_ty == "call":
            return nd1 * self.factor_q
        else:
            return (nd1 - 1) * self.factor_q

    @property
    def gamma(self):
        """ 基础价格对期权价格二阶导 self.ds过小可能不准
        """
        return self.dN(self.d1) * self.factor_q / (self.S * self.sigma * np.sqrt(self.T))

    @property
    def vega(self):
        """波动率对期权价格
        """
        return self.S * np.sqrt(self.T) * self.dN(self.d1) * self.factor_q

    @property
    def theta(self):
        """期限对期权价格
        到期时间不是时间减法,所以要加-
        """
        first = -self.S * self.dN(self.d1) * self.sigma * self.factor_q / (np.sqrt(self.T) * 2)
        if self.opt_ty == "call":
            second = self.q * self.S * self.N(self.d1) * self.factor_q
            third = -self.r * self.K * self.factor_r * self.N(self.d2)
        else:
            second = -self.q * self.S * self.N(-self.d1) * self.factor_q
            third = self.r * self.K * self.factor_r * self.N(-self.d2)

        return first + second + third

    @property
    def rho(self):
        """期限对利率价格
        """
        if self.opt_ty == "call":
            return self.K * self.T * self.factor_r * self.N(self.d2)
        else:
            return -self.K * self.T * self.factor_r * self.N(-self.d2)

    @property
    def vanna(self):
        """ delta对波动率
        """
        return -self.factor_q * self.dN(self.d1) * self.d2 / self.sigma

    @property
    def volga(self):
        """vega对波动率
        """
        return self.factor_q * self.S * np.sqrt(self.T) * self.dN(self.d1) * self.d1 * self.d2 / self.sigma

    @property
    def charm(self):
        """delta对时间,delta衰减率
        """
        first = -q * self.factor_q * self.N(self.d1)
        second = self.factor_q * self.dN(self.d1) / (2 * self.sigma * self.T * np.sqrt(self.T))
        third = 2 * (self.r - self.q) * self.T - (self.d2 * self.sigma * np.sqrt(self.T))
        return second * third + first

    @property
    def speed(self):
        """gamma对现货价格倒数,期权读现货价格三阶导数
        """
        return np.prod([-self.gamma / self.S, 1 + self.d1 / (self.sigma * np.sqrt(self.T))])

    @property
    def color(self):
        """gamma对时间,gamma衰减率
        """
        f1 = np.prod([-self.factor_q, self.dN(self.d1)]) / np.prod([2, self.S, self.sigma, self.T, np.sqrt(self.T)])
        f2 = (2 * (self.r - self.q) * self.T - np.prod([self.d2, self.sigma, np.sqrt(self.T)])) * self.d1 / np.prod(
            [self.sigma, np.sqrt(self.T)])
        f3 = 1 + np.prod([2,
                          self.q,
                          self.T])
        return (f2 + f3) * f1



if __name__ == '__main__':
    # calibrated by https://quantpie.co.uk/oup/oup_bsm_price_greeks.php
    S = 120
    K = 100
    r = 0.03
    opt_ty = "call"
    T = 0.5

    sigma = 0.5
    q = 0.1
    b = BSM(S, K, T, r, opt_ty, sigma, q)
    print("price==", b.price)
    mprice = 50
    print("隐含波动率", b.implied_sigma(mprice))
    b = BSM(S, K, T, r, opt_ty, sigma, q)
    print("delta", b.delta)
    print("gamma", b.gamma)
    print("vega", b.vega)
    print("theta", b.theta)
    print("rho", b.rho)
    print (b.vanna,b.volga,b.charm,b.speed,b.color)

